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GENERATION OF THREE DIMENSIONAL FRACTAL SUBSURFACE 
STRUCTURE BY VORONOI TESSELLATION AND COMPUTATION OF 
GRAVITY RESPONSE OF SUCH FRACTAL STRUCTURE 

Field of the invention 

The present invention relates to a method for the generation of three dimensional 
fractal subsurface structures by Voronoi Tessellation and computation of gravity response of 
such fractal structures. The present invention has a wide range of applications in practically 
all geophysical exploration programs. The method provides the sub-surface basin structure, 
which is very close to the natural setting and hence can be very useful for the modeling of 
hydrocarbon reservoir and mineral deposits for optimum estimation and simulation of the 
reserve having known precisely the extent and shape of the anomalous region by the said 
method, which uses fractal approach. The fractal structure, being very irregular conforms to 
the natural geological situations, which have been utilized herewith to demonstrate its 
applicability in various geophysical exploration programs. 
Background of the invention 

The present invention pertains to provide a new approach for generation of three 
dimensional fractal subsurface structure by Voronoi tessellation and computation of gravity 
response of such fractal structure, which comprises an efficient and entirely new way of 
fractal subsurface generation, which removes all the possibilities of getting into a reentrant 
structure while perturbation of Voronoi centers during iteration steps involved in inverse 
modeling of the underlying structure and computation of gravity response of such fractal 
subsurface, which constitutes forward modeling of the underlying structure for exploration of 
hydrocarbons and minerals. 

It is known from the gravitational law that the bodies having mass would exert 
attractional force on each other. In geophysical aspect instead of mass, the density plays a 
direct role in gravity surveys since density varies laterally as well as vertically, and hence 
affects the mass of a volume. It could be posed in a way that if a large object is having 
density contrast with its surrounding then it will reflect its signature in observed gravity field. 
In case of geophysical studies the density contrast between different interfaces is responsible 
for the gravity anomaly, which in turn can be studied either for hydrocarbon exploration or 
for geological studies also. 

Exploration for hydrocarbons and minerals in subsurface environments with 
anomalous density variations have always been problems for traditional seismic imaging 
techniques by concealing geological structures beneath zones of anomalous density. Many 
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methods for delineating the extent of the highly anomalous density zones exist. Exploration 
for hydrocarbons in subsurface environments with anomalous density variations such as salt 
formations, shale diapers create difficulties in seismic imaging techniques by concealing 
geologic structures beneath zones of anomalous density. By utilizing gravity, magnetic 
measurements and geological constraints along with a robust inversion process, these 
anomalous density zones can be modeled. The spatial resolution obtained from this process is 
normally much lower than that obtained from reflection seismic data. However, models 
obtained from gravity and magnetic data can provide a more accurate starting model for the 
seismic processing. Using the potential field data models as a starting model for two 
dimensional and three dimensional imaging by any technique greatly enhances the 
probability of mapping these concealed geologic structures beneath the zones of anomalous 
density. One source of geologic exploration data that has not been used extensively in spite of 
being cheapest in the past is potential fields data, such as gravity and magnetic data, and 
using potential fields data in combination with seismic data to provide a more accurate depth 
model which is a key parameter in any exploration program. 

In a paper in Geophysical Journal of Royal Astronomical Society, v. 3, 1960, Bott 
suggested a method to trace the floor of a sedimentary basin, which involved the 
approximation of a sedimentary basin by a series of two dimensional juxtaposed 
rectangular/square blocks of uniform density. The assumption of uniform density in entire 
rectangular/ square blocks in said method is highly simplified case of actual reality. 
Computation of gravity anomalies due to two dimensional and three dimensional bodies of 
arbitrary shape has been given in Journal of Geophysical Research, v.64, 1959 by Talwani et 
al. and in Geophysics, v. 25, 1960 by Talwani and Ewing which uses several variables in 
terms of co-ordinates of the corners of the polygons. The method is as accurate as much the 
number of corners is assumed to approximate the irregular shaped body, which enhances 
number of variables to be used in computation. The problem becomes much more 
complicated and cumbersome to deal with so many parameters when the said method is used 
in an inversion algorithm wherein the co-ordinates of the polygon corners are perturbed in 
every iteration in order to fit and achieve the best model. Yet another problem in said method 
arises when used in automated inversion algorithm is that of getting into a reentrant structure 
while perturbation of the polygon corner co-ordinates. This method became very popular and 
is widely used even in present days because of its mathematical simplicity. All these theories 
mentioned above, assumed subsurface structure consisting of some simple geometrical shape 
but Dimri opined the concept of fractal subsurface structure in chapter 16 "Fractal dimension 
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analysis of soil for flow studies" pp. 189-193 In: Application of fractals in earth sciences, 
2000, edited, V. P. Dimri, A. A. Balkema, USA. Benoit Mandelbrot, who coined the term 
fractal in his book "The fractal geometry of nature", W.H. Freeman, New York, 1982 are put 
forth as the fractals are generally irregular (not smooth) in shape, and thus are not objects 
definable by traditional geometry. That means fractals tend to have significant detail, visible 
at any arbitrary scale; when there is self-similarity, this can occur because zooming in simply 
shows similar pictures, alternatively "A geometrical or physical structure having an irregular 
or fragmented shape at all scales of measurement" is known as fractal structure. Earlier in the 
Journal of Geophysical Research v. 100, 1995 Maus and Dimri had opined new concept of 
scaling behavior of potential fields and established a relation between the scaling exponent of 
source and field, which is useful for understanding of fractal geology. The fractal structure 
can be generated using very few parameters by Voronoi tessellation. 

In a straight forward iterative algorithm for the planar Voronoi diagram, Information 
Process Letters, v. 34, 1990, Tipper illustrated that the Voronoi tessellation in a two 
dimensional space consists of enclosing every center by a Voronoi polygon (Fig. 1) such that 
common edge of adjacent polygons is perpendicular bisector to the line joining the centers on 
each side of that edge. Tipper illustrated Voronoi tessellation which uses least square distance 
fromula i.e. L 2 norm which limits the possibility and hence efficiency of generating variety of 
structures merely by changing the exponent p in the L p norm.. Here we have generalised the 
notion of Voronoi tessellation by using L p distances instead of the least square distances as 
mentioned in 'Decon volution and Inverse Theory: Application to Geophysical Problems', 
Elsevier Publication, 1992, by Dimri, V.P., so that Voronoi domains are not necessarily of 
polygonal shape (Fig. 2(b)). 

In yet another study published in proceeding of Indian Academy of Sciences (Earth 
and Planetary Sciences) v. 108, 1999, Moharir et al. have advocated the use of economical 
sub-surface representational techniques for non-local optimization, wherein authors have 
used leminiscate representation which does not provide such an efficient and wider 
possibility of generating subsurface structure as given in the present invention. In the present 
invention a fractal sub-surface is generated which is close to the natural settings and domains 
of different physical property are assigned different colors, which can be later perturbed 
iteratively for inverse modeling, which is important in the interpretation of geophysical data 
as mentioned in the book by Dimri,V.P., 1992, Elsevier Science Publisher, Amsterdam. 

There is a need for a method for efficient and accurate delineation of subsurface 
structure which is close to the real geology present below the earth surface in two dimension 

3 



221/NF/2003 



and three dimension. Such a method should preferably be able to use physical property 
variations in the subsurface geometry. Also, the method should preferably be able to obtain 
realistic forward model of anomalous formations that are expected in normal sedimentary 
basin, which can be used in exploration. The present invention satisfies this need. 
Objects of the invention 

The main object of the invention is to provide a method for the generation of three 
dimensional fractal subsurface structures by Voronoi Tessellation and computation of gravity 
response of such fractal structure, which obviates the drawbacks as detailed above. 

Another object of the present invention is to provide an efficient method for 
generation of fractal subsurface which is very close to real geological situations. 

Yet another object of the present invention is to show that the present invention has its 
direct implication in reservoir modeling by way of generating complex geological situations 
wherein the variation of physical properties can be very well studied and predicted using 
fractal approach. 

Yet another object of the present invention is to provide an excellent idea which can 
be used any branch of geophysical exploration. 
Summary of the invention 

Accordingly the present invention provides a method for the generation of three 
dimensional fractal subsurface structure by Voronoi tessellation and computation of gravity 
response of such fractal structure, which comprises a robust and efficient process for 
generation of fractal subsurface structures, which is very close to natural setting of the 
subsurface geology and provides computation of forward gravity response of such structures 
for delineation of the underlying anomalous object. 

In an embodiment of the present invention Voronoi centers may be selected and the 
region of interest may be tessellated around those centers to provide a geologic subsurface 
earth model with physical property variation viz. density in case of gravity prospecting. 

In another embodiment of the present invention the gravity response due to said 
geological subsurface earth model is computed, which can be introduced into the inversion 
process to further refine the model or as is done in forward modeling, the model may be 
changed by perturbing Voronoi centers and again gravity response may be computed to 
match the observed gravity anomaly. This process of inversion /forward modeling is repeated 
until the results converge to a single answer. 

In yet another embodiment of the present invention the method is used to provide an 
excellent technique which has its application in various geophysical methods other than 
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gravity e.g. seismic for the precise mapping of reservoir geometry which intern shall help in 
demarcating reservoir boundary, physical property variations and as a whole in enhancing oil 
recovery from an existing oil reservoir. 
Brief description of the accompanying drawings 

Fig. 1 (a) and (b) represents the fundamental of existing Voronoi tessellation where in 
circles are drawn (shown in blue color) with increasing radius about the Voronoi centers 
(shown in red color) and a tangent line (shown in green color) is drawn where the adjacent 
circles touches each other, (b) shows the particular case of generating layered model. 

Fig. 2 (a) Represents a general case where a polygonal structure is generated using 15 
Voronoi centers (shown as red dots) using existing Voronoi tessellation, (b) represents a 
subsurface generated by modified Voronoi tessellation using four Voronoi centers (shown as 
black dots). 

Fig. 3 represents the 3-D subsurface structure wherein 2-D layers of variable physical 
property regions are overlain to generate 3-D volume. 2-D layers are shown separate for the 
sake of clarity of pyhsical property variations shown in each layer. The figure demonstrate 
the ability of generating various kind of structures shown in different layers merely by 
changing exponent p in L p norm, keeping Voronoi centers fixed. 

Fig.4. represents the 3-D structure wherein layers of different physical property 
variations are shown. The figure demonstrate the another possibility of generating various 
structures as shown in different layers by changing Voronoi centers keeping exponent p 
constant in L p norm. 

Fig. 5 represents grid overlain on the region of interest. 

Fig. 6 (a) represents an example of gravity anomaly response over the simplified 
horizontal layered model, (b) Represents an example of gravity response calculated over 
fractal subsurface structure showing four regions of different physical property (density) 
variations. Density values assigned to red, blue, green and magenta colors are 2.1, 2.3, 2.67 
and 2.5 g/cc respectively, and the depth at which fractal subsurface is assumed is 10 units. 

Fig. 7 represents flowchart of the procedure adopted for the computation of gravity 
response over fractal subsurface. 
Detailed description of the invention 

The present invention brings a new facet of domain characterization by a set of 
parameters, referred herein as Voronoi centers. These parameters are perturbed and thus the 
different tessellated regions are generated at different depths. This characterization is very 
suitable and entails the development for solution of geophysical inverse problems with the 
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help of non-local optimization techniques. 

Given a set of points qi,q 2 , q n in R n (Fig. 2 (a), points are shown as red dots), the Voronoi 

Region Vj of a particular vector q { is defined as the set of all points in R n for which q* is the 
nearest vector: 



where R n is set of real numbers, q* are the points inside the Voronoi region V i? x is any point 

within Vj <z R n , N is any integer. 

In a straight forward iterative algorithm for the planar Voronoi diagram, Information 

Process Letters, v. 34, 1990, Tipper illustrated that the Voronoi tessellation in a two 

dimensional space consists of enclosing every center by a Voronoi polygon (Fig. 1) such that 

common edge of adjacent polygons is perpendicular bisector to the line joining the centers on 

each side of that edge. Here we have generalised the notion of Voronoi tessellation by using 
p 

L distances instead of the least square distances so that Voronoi domains are not necessarily 
of polygonal shape (Fig. 2(b)). The L p distance is given by the expression: 

If =(x — q j ) 1/p , where x is an arbitrary point and qj is a point whose distance has to be 

calculated, and p is an exponent which can hold any real value. Initially a few Voronoi 
centers are taken and from them two dimensional Voronoi tessellation is generated in which 
different domains are shown by different colours (Fig. 2(b)). Density contrast corresponding 
to each domains of the tessellated region is assigned while its generation by Voronoi 
tessellation. A grid is laid over the tessellated region on the surface and gravity anomaly is 
calculated at each grid node (Fig. 5). The gravity response due to each subsurface region of 
different physical property (shown by different colors in Fig. 2(b)) is calculated at each node 
of the grid laid at the surface. Next the cumulative response of each subsurface region of 
different layers (Fig. 3) is computed at each nodes of the grid laid at the surface. If tessellated 
regions are at equal depths then the integrated response can be calculated by Simpson's rule 
other wise Gauss's quadrature formula can be used for tessellated regions at unequal depths. 
The gravity anomaly at each nodes of the grid laid at surface is computed due to subsurface 
tessellated region at depth Z=10 unit as an example shown in Fig. 6. 
Mathematical expressions for the computation of gravity anomaly 

The gravity anomaly caused by the polygonal lamina per unit thickness, in a form 
suitable for programming is expressed in terms of Xi,y i5 z and x i+ i,y i+ i,z, the co-ordinates of 
two successive vertices of the polygon and has been given in Geophysics, v. 25, 1960, by 
Talwani and Ewing as: 
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V = GpUW arccosfc I r,)(x l+1 I r M ) + (y i /r i )(y i+l I r M )} 



- arc sin — z — * . + arc sin 



1/2 



/■ 2 . 2x1/2 / 2 , 2x 

(Pi + Z, ) (Pi + Zi ) 



Where S= +1 if pi is positive, S= -1 if pi is negative, 

W= +1 if mi is positive, W= -1 if mj is negative, 

'Z' is depth and 'n' is number of sides in the polygon. 

G is universal gravitational constant, p is the density of the tessellated regions (shown in 
diffeerent colours in Fig. 2(b)). 

Pi = — x t — l - ^ y t 9 

Hi ' j 

y = x / ~~ ^+1 jg+i + >j ~ yj+i y e -+i 

r i+\ r i 9 i+l r i+\ 

r i+l 'j r i+l r » 
, / 2 . 2x1/2 

r i+l = + (x i+l 2 + y i+1 2 ) 1 ' 2 , 

r i , i+l =H(x i -x i+l ) 2 +(y i -y i+l ) 2 ] u2 . 

In the application of this algorithm the effect of common arm adjacent to the 
polygonal tessellated region is removed as it contributes twice for the response in gravity 
anomaly computation. The stepwise implementation of the algorithm is shown in Fig. 7 as a 
flowchart. 

The present invention has its novelty over previous work in following counts: 

• The method uses fractal subsurface structure which represents complex geological 
structures. Simple polygonal structures were used because of mathematical simplicity 
and such realistic structures were avoided hitherto. 
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• The method uses physical property variation in lateral and vertical direction; provides 
the solution more closure to reality because none of the things in general are found 
homogenous beneath the earth as has been assumed in previous literature. 

• The method is very efficient as it involves very few Voronoi centers to change the 
geometry or even keeping the Voronoi centers constant we are able to achieve various 
combinations of the complex structures merely by changing the exponent p in L p 
norm used in Voronoi tessellation (Fig. 4). This is in contrast to the previous work 
where simple polygonal geometry was used and to represent polygons one has to deal 
with as many parameters as the number of corners in the polygons. 

• The present method is useful for delineation of complex subsurface structures, which 
are generally favourable for oil accumulation. The previous work done in this line 
may not be as accurate and efficient as the present work because of the simple 
geometry used for the modeling of the underlying objects. 

The following examples are given by way of illustration and therefore should not be 
construed to limit the scope of the present invention. 
EXAMPLE 1 

A fractal subsurface is generated using L p norm taking p=1.5, wherein the Voronoi region 
is defined by the co-ordinates x= 20 to 45 and y= 25 to 50. The Voronoi centers chosen 
within this region are given by the following co-ordinates: 



X 


y 


22.6 


28.75 


30.6 


30.6 


36.6 


36.6 


41.4 


41.4 



The subsurface thus generated is shown in Fig. 2(b), with Voronoi centers marked as 
black dots. 
EXAMPLE 2 

In yet another example the fractal subsurface at different depth levels are generated 
using different L p norms. The co-ordinates of Vornoi region for all the subsurfaces were 
taken as x= 20 to 50 and y= 5 to 65 and Voronoi centers within the region were taken as: 



X 


y 


22.0 


26.0 


30.0 


33.0 1 


37.0 


40.0 


46.0 


48.0 


48.0 


58.0 
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This example shows the variation in geometry and provides and excellent way of 
changing the geometry merely by changing the exponent p in L p norm. The results are shown 
in Fig. 3 where the first layer (topmost) corresponds to p=1.5, the middle layer corresponds to 
p=-1.5 and the lowermost corresponds to p=1.0, which is equivalent to L 1 norm. 
EXAMPLE 3 

In yet another example as shown in Fig. 4 the fractal subsurface at different depth 
levels are generated using same L p norm, where p=1.5 but Voronoi centers are changed. This 
example illustrates the possibilities of generating different kind of structures by changing 
Voronoi centers. The co-ordinates of Vornoi region for all the subsurfaces were taken as x= 
20 to 50 and y= 5 to 65 and 5 Voronoi centers within the Voronoi region were taken having 
following x,y co-ordinates: 



Top layer 


Middle Layer 


Bottom Layer 


xl 


yi 


x2 


y2 


x3 


y3 


22.0 


26.0 


22.0 


7.0 


15.0 


8.0 


30.0 


33.0 


22.0 


50.0 


30.0 


50.0 


37.0 


40.0 


35.0 


33.0 


35.0 


10.0 


46.0 


48.0 


50.0 


7.0 


45.0 


45.0 


48.0 


58.0 


50.0 


50.0 


55.0 


10.0 



EXAMPLE 4 

In yet another example the gravity anomaly computed over the fractal subsurface is 
shown in Fig. 6. Figure 6(a) corresponds to the particular case wherein Voronoi centers are 
taken along a line which gives layered model where as Fig. 6(b) represents the fractal 
structure where Voronoi centers correspond to those given in EXAMPLE 1. The density 
values corresponding to red, blue, green and magenta colors are 2.1, 2.3, 2.67 and 2.5 
respectively. The fractal subsurface is assumed at the depth of 10 unit from the surface. 
THE MAIN ADVANTAGES OF THE INVENTION ARE: 

1. To provide a new approach for generation of three dimensional fractal subsurface 
structure by Voronoi tessellation and computation of gravity response of such fractal 
structure, this is an efficient way of fractal subsurface generation and computation of 
gravity response of the same. This method removes all the possibilities of getting into a 
reentrant structure while perturbation of Voronoi centers during iteration steps involved 
in inverse modelling of the underlying structure. The said method constitutes delineation 
of the underlying structure for exploration and production of hydrocarbons and minerals 
lying below the earth's surface. 
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2. The technique of Voronoi tessellation has been modified using L p norm, where in p can 
assume any value from the set of rational numbers, provides wider possibilities of 
generating different realistic sub-surface structures. This is an innovative approach 
proposed herewith for delineation of complicated geologic structure which are favorable 
for oil accumulation and was never attempted by any one hitherto. 

3. The method provides the sub-surface basin structure, which is very close to natural setting 

and hence can be very useful for the modeling of hydrocarbon reservoir and mineral 
deposits for optimum estimation and extraction of the reserve having known precisely the 
extent and shape of the anomalous region by the said method. This remained un- 
attempted in the earlier work. 

4. The modified Voronoi tessellation as described herewith generates fractal subsurface 
merely by changing Voronoi centers in the inverse modeling process, unlike perturbation 
of co-ordinates of all the vertices of the geometrical structure used in all the existing 
modeling methods. This makes the process very efficient and fast and also provides wider 
possibilities of efficient modeling. 

5. A new approach for generation of three dimensional fractal subsurface structure by 
Voronoi tessellation and computation of gravity response of such fractal structure is 
entirely a new invention which has got great implication in discovery of hydrocarbon and 
minerals. 
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